function [ F ] = F_distribution_power(q, q_bar, qmin, betaq, sig, q_up)

if false
    q = [8:0.1:14];
    qmin = 0.01;
    betaq = 1.15;
    sig = 0.01;
    q_up = 13;
end

c = (betaq-1)/(qmin^(1-betaq));

if sig == 0
    F = (q < q_bar) .* (c./(1-betaq)).* (q.^(1-betaq)- qmin^(1-betaq))+ ...
        (q>=q_bar) .* (q<q_up) .* (c./(1-betaq)).* (q_up.^(1-betaq)-qmin^(1-betaq)) + ...
        (q>=q_up) .* (c./(1-betaq)).* (q.^(1-betaq)-qmin^(1-betaq));
else
    F = zeros(length(q),1);
     
    for i = 1:length(q)
 
        F(i) = integral(@(u) f_density_power(u,q_bar,qmin,betaq,sig,q_up),qmin,q(i), 'ArrayValued',true);
        
    end
end
end

